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Abstract 

Response surface functions are often used as simple and inexpensive replacements for computationally 
expensive computer models that simulate the behavior of a complex system over some parameter space. 
“Progressive” response surfaces are ones that are built up progressively as global information is added from 
new sample points in the parameter space. As the response surfaces are globally upgraded based on new 
information, heuristic indications of the convergence of the response surface approximation to the exact (fit- 
ted) function can be inferred. Sampling points can be incrementally added in a structured fashion, or in an 
unstructured fashion. Whatever the approach, at least in early stages of sampling it is usually desirable to 
sample the entire parameter space uniformly. At later stages of sampling, depending on the nature of the 
quantity being resolved, it may be desirable to continue sampling uniformly over the entire parameter space 
(Progressive response surfaces), or to switch to a focusing/economizing strategy of preferentially sampling 
certain regions of the parameter space based on information gained in early stages of sampling (Adaptive 
response surfaces). Here we consider Progressive response surfaces where a balanced indication of global 
response over the parameter space is desired. We use a variant of Moving Least Squares to fit and interpolate 
structured and unstructured point sets over the parameter space. On a 2-D test problem we compare response 
surface accuracy for three incremental sampling methods: Progressive Lattice Sampling; Simple-Random 
Monte Carlo; and Halton Quasi-Monte-Carlo sequences. We are ultimately after a system for constructing 
efficiently upgradable response surface approximations with reliable error estimates. 

Introduction and Background 

Large-scale optimization and uncertainty analyses are often made feasible through the use 
of response surfaces as surrogates for computational models that may not be directly 
employable because of prohibitive expense and/or noise properties and/or coupling diffi- 
culties in multidisciplinary analysis. Examples of response surface usage to facilitate 
large-scale optimization and uncertainty analyses are cited in Roux et al. (1996), Unal el 
al. (1996), and Venter et al. (1996). 

Two issues that arise when using response surface approximations (RSA) are accuracy and 
the cost of procuring the data samples needed to create the RSA. With a sufficiently flexi- 
ble global fitting/interpolating function over the parameter space, response surface accu- 
racy generally increases as the number of data points increases (if the points are 
appropriately placed throughout the parameter space), until the essential character of the 
function is effectively mapped out. Thereafter, it is not cost effective to continue adding 
samples. Since a single high-fidelity physics simulation (i.e., one data sample) can take 
many hours to compute, it is desirable to minimize the number of simulations that are 
needed to construct an accurate response surface. 

For our purposes here it is assumed that: 1) the computer model is relatively expensive to 
evaluate; 2) the parameter space is a unit hypercube or can be accurately and inexpen- 
sively mapped into one; 3) the sampled or “target” function is a continuous, deterministic 
function over the parameter space; 4) reasonably general, arbitrary target functions are to 
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be fitted; and 5) approximate response values are desired over the entire parameter space 
or subspace being considered -i.e., for global and local optimization or mapping inputs to 
outputs in uncertainty propagation. 

Given these specifications, Romero et at. (2000) examined several formulations for con- 
structing progressive response surfaces built on Progressive Lattice Sampling (PLS) incre- 
mental sampling designs. PLS is a paradigm for structured uniform sampling of a 
hypercube parameter space by placing and incrementally adding sets of samples such that 
all samples are efficiently leveraged as the design progresses from one level to the next. 
Figures 1-4 show successive PLS levels in 2-dimensions. (Also shown for comparison 
are point sets from classical simple-random Monte Carlo sampling -using three different 
seeds for the random-number generator, and from Halton “quasi- Monte Carlo” low-dis- 
crepancy sequences (see, e.g., Owen, 2003) where we have used two different sets of 
prime-number bases to create the two sets of Halton samples.) 

PLS endeavors to preserve uniformity of sampling coverage over the parameter space in 
the various stages or levels of the incremental experimental design. Uniform coverage 
over the parameter space is desirable for general response surface construction because 
this reduces the redundancy or marginalization of new information from added samples. 
This is a basic concept of upgradable quadrature methods (Patterson- 1968, and Genz & 
Malik- 1983). 

PLS builds knowledge by reducing global knowledge deficit over the parameter space. It 
does not attempt to build specific or targeted knowledge by building on previous informa- 
tion in the manner of “adaptive” sampling, which efficiently maximizes knowledge over 
particular regions of the parameter space. Thus, PLS designs select sample locations 
strictly on geometric principles such that new samples are intended to be “maximally far” 
from each other and from all other existing samples at each level of sampling. Thus, global 
uniformity of coverage is maintained at each level as the scheme progresses. 

The arrangement of samples in each PLS level allows the parameter space to be subdi- 
vided into a regular pattern of adjacent polygons, which for two parameter space dimen- 
sions results in triangular and quadrilateral 2-D finite elements (FEs) yielding linear to 
quadratic polynomial interpolation over each element (see Romero & Bankston, 1998). 
The collection of all the elements together creates a C°-continuous global function over 
the parameter space. As such, the global RSA has considerable freedom to locally con- 
form to the data values of the sample points (see Figures 6 - 8). A mathematical analysis 
(Strang & Fix, 1973) of finite element interpolation of this nature shows that for a continu- 
ous and infinitely differentiable function over the parameter space, the domain integral of 
the pointwise absolute error goes to zero as the spacing between samples goes to zero. 
This analysis can be applied to the FE/PLS method. Hence, in the limit of infinite sam- 
pling, FE/PLS response surfaces converge everywhere in the domain to target functions in 
this class (though the convergence rate is generally not uniform over the domain). There- 
fore the FE/PLS method provides a convergent reference against which the accuracy of 
other progressive response surface methods can be compared. 

Upon upgrading from one level to the next of the PLS design, the resulting change in the 
response surface at any point in the domain is a heuristic indicator of the magnitude of the 
local error in the response surface approximation. When the incremental change goes to 
zero, this tentatively indicates that the local approximation error has become negligible. 
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SRS9-seedl SRS9-seed2 SRS9-seed3 



Halton9-2,3 Halton9-5,7 PLS-Level4 



Figure 1. 9-point sample sets on a 2-D unit Hypercube from three types of incremental 
sampling methods: A) Top Row- SRS Monte Carlo with three different initial 
seeds; and B) Bottom Row- base 2,3 and 5,7 Halton sets, and PLS set. 


SRS13-seedl SRS13-seed2 SRS13-seed3 



Haltonl3-2,3 Haltonl3-5,7 PLS-Level5 



Figure 2. 13-point sample sets on a 2-D unit Hcube from three types of incremental sam- 
pling methods: A) Top Row- SRS Monte Carlo with three different initial seeds; 
and B) Bottom Row- base 2,3 and 5,7 Halton sets, and PLS set. 
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Figure 3. 25-point sample sets on a 2-D unit Hcube from three types of incremental sam- 
pling methods: A) Top Row- SRS Monte Carlo with three different initial seeds; 
and B) Bottom Row- base 2,3 and 5,7 Halton sets, and PLS set. 
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Figure 4. 41-point sample sets on a 2-D unit Hcube from three types of incremental sam- 
pling methods: A) Top Row- SRS Monte Carlo with three different initial seeds; 
and B) Bottom Row- base 2,3 and 5,7 Halton sets, and PLS set. 
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The advantage of the structured PLS approach is that it is thought to be globally optimal or 
nearly optimal in that, as samples are added in attaining each new Lattice Sampling 
“Level”, the spacing of samples throughout the parameter space remains uniform, or 
nearly so, given the constraint of previous sample locations. The locations of previous 
samples are respected because it is desired to fully leverage them (with minimal marginal- 
ization) as new samples are added. 

This would appear to be the most efficient way (following 
the precedent of upgradable quadrature methods) to pro- 
gressively build up a response surface. It may be that at any 
given level if complete freedom is allowed where the sam- 
ples can be placed, then these may be arrangeable over the 
parameter space with better uniformity than PLS provides. 

For example, Romero et al. (2003) find that this appears to 
occur sometimes with Latin Hypercube Sampling and Cen- 
troidal Voronoi Tessellation -depending on the random 
number generator initial seed, number of samples, etc. 

However, these are non-incremental sampling methods; 
augmenting the number of samples would imply a com- 
pletely different sampling of the parameter space at all new 
point locations. To go from X to X+Y samples in the space 
would therefore involve X+Y new evaluations of the gener- 
ating function. This is more expensive than an incremental 
sampling method like PLS, which at each new stage costs 
only the increment of Y new samples (to be added to the original X for a total of X+Y uni- 
formly dispersed samples). 

A strong disadvantage of PLS, however, is that its structured experimental-design sam- 
pling nature allows only a quantized increment Y to be added to an existing PLS level 
(point set) to graduate to a new level. Hence, there is a constraint on the number of sam- 
ples that can be added and still maintaining as uniform a filling of the hypercube as the 



Figure 5. 2-D test func- 
tion to be 
approximated 
(“generating” 
function). 



(A) (B) (C) 


Figure 6. Successive response surface approximations to the exact generating func- 
tion by finite-element interpolation of PLS designs: (A) - Level 1 having 3 
samples; (B) - Level 3(4) having 9 samples; (C) - Level 6 having 25 samples. 
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scheme is capable of. In particular, this quantized incremental cost Y scales increasingly 
quickly as the PLS level and dimension of the parameter space increase. 

Other, incremental sampling methods exist that are unstructured (non-quantized) that 
don’t have the hard cost-scaling of PLS. Some popular unstructured incremental- sampling 
methods are the Halton “quasi- Monte Carlo” low-discrepancy sequence method (Owen, 
2003), and Simple-Random Sampling (SRS) Monte Carlo. These allow Y as small as 1 
without creating a prejudiced imbalance in the global coverage/filling of the parameter 
space. Hence, regardless of the dimension of the space and the stage of sampling, we 
could incrementally add, say, 50% more samples at a time in monitoring the convergence 
of the response surface results. If the sampling budget is reached before adequate conver- 
gence is established, then the final results will still be based on a point placement that is 
characteristic of the most uniform sample placement the method is capable of producing. 

The point placement of Halton sampling appears characteristically more uniform over the 
parameter space than that of SRS, as can be seen in Figures 1 - 4. (We assert this based 
only on graphical appearance, not on formal measures of uniformity. Burkhardt et al. 
(2004) describe reliable metrics of point uniformity in hypercubes, but these have not been 
applied here.) In turn, the uniformity of PLS appears to be better than that of Halton. In the 
next section we will confirm that better sampling uniformity over the parameter space gen- 
erally correlates with better response surface accuracy. 

Given a set of sampling points over a parameter space, the quality of the response surface 
approximation (RSA) also depends on the particular method used to fit and interpolate the 
data. We now turn to consideration of data fitting methods to interpolate and extrapolate 
the sample data. Finite-element interpolation of sample data for general unstructured point 
placement and arbitrary numbers of samples in arbitrary dimensions is a difficult prospect. 
It is not immediately obvious how anything but linear tetrahedral (simplex) elements could 
readily be used, thereby sacrificing any higher-order convergence potential in the piece- 
wise-linear interpolation scheme. Extrapolation procedures are also not immediately obvi- 
ous, as this is not normally encountered in the Finite Element Method. 

Four general data fitting and interpolation/extrapolation methods that can work with struc- 
tured or unstructured progressive sampling schemes have been evaluated by the authors. 
These are global polynomial regression and kriging (Romero, et al., 2000, and Krishna- 
murthy & Romero, 2002); moving least squares (Krishnamurthy & Romero, 2002, and 
Romero et al., 2003); and Radial Basis Function methods (Krishnamurthy, 2003). Though 
our experience is very limited, of these, MLS appears to present a good balance of 
response surface accuracy, smoothness, robustness, and ease of use. Therefore, we use 
MLS in the next section to generate response surfaces from the 2-D sample sets in Figures 

1- 4. 

2- D Example Problem for Examining Performance of Progressive Response 
Surfaces 


Figure 5 plots a 2-D model function used to study the effect of sample point addition, 
sample placement scheme, and interpolation method on response surface accuracy. The 
value of the multimodal function in Figure 5 is defined as: 


Y( pl,p2)= 0.8r + 0.35 sin 


2.4:i-L 


72 


[1.5 sin (1.30)] 


(EQ 1) 
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where r = J(pi ) 2 + (p2 ) 2 and 0 = arctan^j on the domain o<pi<i,o<p2<i. 

Exact values (samples) of this function are obtained at 9, 13, 25, and 41 points shown in 
Figures 1 - 4 for the PLS, SRS, and Halton incremental sampling methods. These numbers 
of samples correspond to PLS levels 4 - 7 in 2-D. (The high deterministic uniformity of 
the PLS sets provides a good standard to compare the other point sets against.) 

The point data is then fitted and interpolated with Moving Least Squares (MLS). The par- 
ticular implementation of MLS we use is described in Krishnamurthy & Romero (2002). 
The quintic weight function described in the reference is used here to give C2 smoothe- 
ness to the MLS global interpolation function over the 2-D parameter space. A quadratic 
polynomial basis function is used for local interpolation, which requires at least 
(M+l)(M+2)/2 sample points (6 for M=2 dimensions) within a given evaluation point’s 
local radius of influence. An optimal local radius of influence is calculated and used for 
each different point set, so that this element of fitting error is minimized in this study. 

In order to assess the response surface error due to the fitting method (as opposed to 
the number and location of sample points), the finite element interpolation method of 
Romero & Bangston (1998) is also used to fit the PLS data sets. Figure 6 shows FE/PLS 
response surfaces for sets of 3, 9, and 25 samples. This illustrates the convergence of the 
progressive response surface to the target function as samples are added from the PLS 
design. 

To examine the fitting performance of the progressive response surfaces, a simple 
measure of quality of global fit over the parameter space is used: 

441 

^ [exact- - predicted^ 

approximate spatial average absolute error = — (EQ 2) 

where exact and predicted values in the summation come from respective evaluation of the 
exact function and the particular response surface approximation at 441 equally spaced 
points on a 21x21 square grid overlaid on the unit-square domain. This measure is an 
expedient approximation to the global average integrated absolute error over the domain, 
which would require a much more involved calculation. 

Performance of Progressive Response Surface Methods on 2-D Problem 

Figure 7 presents the response surface errors for the various sampling schemes, numbers 
of samples, and interpolation methods. All sampling methods yield a reduction of global 
fitting error as the number of samples increases in progressing to each new stage of the 
response surface. 

At every population level (9, 13, 25, 41 samples), SRS-based sample placement performs 
worst (for all three of the random-number- generator initial seeds tried). This is a general 
reflection of the less uniformity with which the points are placed in the domain, as Figures 
1-4 reveal. The more uniform Halton placement (for both the 2,3 and 5,7 prime bases for 
pl,p2 coordinates) performs significantly better than SRS in general, but still significantly 
less well than the deterministically uniform PLS point placement. The effect of point 
placement is substantial; at a population level the error difference between the best and 
worst point placements is roughly equal to the error difference between that and the next 
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population level. That is, using a better point-placement scheme can substantially reduce 
the number of samples (and evaluation cost) needed to achieve a given level of response 
surface fitting error. 


Among the incremental sampling methods (which are the most cost efficient type of sam- 
plers for progressive response surface construction), of the ones tried here Halton appears 
to be the best choice in general. Being one-at-a-time incremental sampling methods, Hal- 
ton and SRS don’t have the hard cost-scaling of PLS (where only quantized increments are 
allowed), and Halton point uniformity is much better than SRS uniformity in general. 


With regard to interpolation methods, both FE and MLS are applied to the PLS sets for 
comparison. Which interpolator is used has relatively less effect here than the number and 
placement of sample points, but the effect is still non-negligible at 25 and 41 points. 
Though FE interpolation generally yields smaller error than MLS, this is not always true 
(e.g., the data for 13 samples). Futhermore, FE is not currently a viable choice for unstruc- 
tured Halton sampling, which is the most viable of the incremental sampling methods con- 
sidered here. 


Sampling types vs. average error 
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Figure 7. Convergence of progressive response surfaces to generating function 
as various sampling methods add samples. (Spatially averaged absolute 
error is plotted.) 


Closing 

It can be very difficult to determine when a particular sampling design and interpola- 
tion scheme sufficiently resolve a function, yet this must be done if the response surface is 
to be used as an effective inexpensive replacement for the actual function. Monitoring con- 
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vergence heuristics of progressive convergent response surface approximations can help in 
this regard, and the Halton/MLS combination initially appears to be a generally viable 
approach for generating such response surfaces. However, there are still some subtle 
issues with both of these technologies that have to be characterized and addressed before 
they can be ready for general robust implementation. This will be the subject of future 
development, testing, and evaluation of the Halton/MLS method (and any other viable 
alternatives identified in the future). 
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